Consider the simple-looking second-order differential equation

$\frac{d^{2}g}{d{t}^2}+g^n=1$

where the exponent n can take any positive real value. This equation appropriately scaled describes the evolution of the square of the radial variable of a body moving under the influence of a power potential.

The equation is in general hopelessly nonlinear and as such not amenable to exact solution. It turns out, however, that a very good approximate solution can be found for continuous variation of the exponent. Initial conditions for the solution will have a vanishing first derivative at the origin, which restricts the solution to functions even in the variable of evolution. This corresponds to a body beginning movement from one turning point.

Numerical solutions with these initial conditions look very similar as a function of the exponent, as can be readily seen in this interactive graphic:

The initial condition c₀ is limited by the value n of the exponent in a manner to be described later.

For all values of the exponent there is one trivial exact solution: $g\left(t\right)=1$ . This solution corresponds to a circular orbit at the maximum angular momentum allowed for bound orbits.

For positive values of the exponent there are two simple exact solutions. For n = 1 the solution is a simple sinusoidal,

$g\left(t\right)=1+(c_0-1)cost,n=1$

as can be verified by differentiation. The solution is symmetric about the trivial solution, with extrema at $1+(c_0-1)$ and $1-(c_0-1)$ . Using a half-angle identity this can also be written

$g\left(t\right)=2-c_0+2(c_0-1){cos}^2\frac{t}{2},n=1$

For n = 2 the exact solution is given by elliptic functions. Since any elliptic function can be expressed in terms of other elliptic functions with the same periods, choose a representation of the solution that relates easily to the form just given:

$g\left(t\right)=c_1+(c_0-c_1){cn}^2(\omega t,k),n=2$

The solution is no longer symmetric about the trivial solution, necessitating two constants for the extrema. The second of the two extrema occurs at quarter periods of the Jacobi elliptic function,

$\omega t_1=K\left(k\right)\to cn\left[K\right(k),k]=0$

where the elliptic function becomes zero. The second derivative of the squared Jacobi elliptic function is

$\begin{array}{l}\frac{d^2}{d{t}^2}{cn}^2(\omega t,k)=2{\omega }^2[{sn}^2(\omega t,k){dn}^2(\omega t,k)-{cn}^2(\omega t,k){dn}^2(\omega t,k)\\ \hfill +k^2{sn}^2(\omega t,k){cn}^2(\omega t,k)]\\ \phantom{\frac{d^2}{d{t}^2}{cn}^2(\omega t,k)}=2{\omega }^2[1-k^2+2(2k^2-1){cn}^2(\omega t,k)-3k^2{cn}^4(\omega t,k)]\end{array}$

Substituting this in the differential equation leads to a set of equations among the four constants,

$\begin{array}{r}2{\omega }^2(1-k^2)(c_0-c_1)+c_1^2=1\\ 2{\omega }^2(2k^2-1)+c_1=0\\ -6{\omega }^2{k}^2+(c_0-c_1)=0\end{array}$

Using the last two equations to eliminate the extrema from the first equation gives

$4{\omega }^4(k^4-k^2+1)=1$

while adding the last two equations together gives the relation

$2{\omega }^2(k^2+1)=c_0$

Eliminating the frequency between these two relations, the elliptic modulus is a solution of

$\begin{array}{c}(c_0^2-1)k^4-(c_0^2+2)k^2+c_0^2-1=0\\ k^2=\frac{c_0^2+2\pm \sqrt{3}{c}_0\sqrt{4-c_0^2}}{2(c_0^2-1)}=\frac{[\sqrt{3}{c}_0\pm \sqrt{4-c_0^2}{]}^2}{4(c_0^2-1)}\end{array}$

The square of the frequency is

${\omega }^2=\frac{c_0}{2(k^2+1)}=\frac{c_0^2-1}{\sqrt{3}[\sqrt{3}{c}_0\pm \sqrt{4-c_0^2}]}$

and the second extremum is

$c_1=c_0-6{\omega }^2{k}^2=c_0-\frac{\sqrt{3}}{2}[\sqrt{3}{c}_0\pm \sqrt{4-c_0^2}]$

Since both extrema are equal to unity for the trivial exact solution, the correct sign for the radical is the lower one. The final solution, with a bit of simplification, is thus

$\begin{array}{l}{c}_1=-\frac{c_0}{2}+\frac{3}{2}\sqrt{\frac{4-c_0^2}{3}}\\ \omega =\frac{1}{2}\sqrt{c_0+\sqrt{\frac{4-c_0^2}{3}}}\\ k=\frac{\sqrt{3}{c}_0-\sqrt{4-c_0^2}}{2\sqrt{c_0^2-1}}=\frac{2\sqrt{c_0^2-1}}{\sqrt{3}{c}_0+\sqrt{4-c_0^2}}\end{array}$

The frequency and elliptic modulus have been rationalized to show that they both remain finite for the trivial solution, with the latter going to zero.

For an solution analytic in the exponent n, the constants here would include a functional dependence on that value. Since it is a property of the Jacobi elliptic function that $cn(x,k)\to cosx$ for zero modulus, the exact solution for n = 2 has a simple relationship to that for n = 1 . If the functional dependence on exponent of the former solution were to allow k → 0 for arbitrary c₀ as n → 1 , then the exact solution for n = 1 would then be recovered with the additional limits

$c_1\to 2-c_0\mathrm{and}\omega \to \frac{1}{2},n\to 1$

While such a solution analytic in the exponent n is not forthcoming, the relationship between these two exact solutions suggests that one of them can be used as an approximate solution for all values of the exponent. Since the the squared Jacobi function already displays an asymmetry about the trivial solution that the circular function does not, the former is a more appropriate choice for the approximation.

Of the three constants needed apart from the initial condition to describe the approximate solution, two can be stated exactly. First multiply the original differential equation by the first derivative of the dependent variable and integrate once:

$\begin{array}{c}\frac{dg}{dt}[\frac{d^{2}g}{d{t}^2}+g^n=1]\\ \frac{1}{2}(\frac{dg}{dt}{)}^2+\frac{g^{n+1}}{n+1}=g-l\end{array}$

The constant of integration l is proportional to the square of angular momentum in physical applications, hence the letter. At the extrema of the motion the first derivative goes to zero so that

$c_0-\frac{c_0^{n+1}}{n+1}=l=c_1-\frac{c_1^{n+1}}{n+1}$

The second extremum is thus a solution to the non-polynomial equation

$\frac{c_1^{n+1}}{n+1}-c_1+l=0$

With the scaling $c_1=(n+1{)}^{1/n}c$ this becomes

$c^{n+1}-c+\frac{l}{(n+1{)}^{1/n}}=0$

which is close to the functional form given at the beginning of a presentation on the real roots of non-polynomials. Accounting for minor differences, the solution found for the root in the vicinity of zero is

$c_{\mathrm{zero}}=\sum _{i=0}^{\infty }\frac{\Gamma [ni+i+1]}{\Gamma [ni+2]}\frac{1}{i!}[\frac{l}{(n+1{)}^{1/n}}{]}^{ni+1}$

and the second extremum in terms of the first is

$c_1=(n+1{)}^{1/n}\sum _{i=0}^{\infty }\frac{\Gamma [ni+i+1]}{\Gamma [ni+2]}\frac{1}{i!}[\frac{(n+1)c_0-c_0^{n+1}}{(n+1{)}^{1+1/n}}{]}^{ni+1}$

It was demonstrated in the presentation that this series converges for its entire domain of application, so that in principle it can be taken as an exact result. In practice it is simpler and more accurate to evaluate the second extremum as a numerical root of the non-polynomial equation.

The real periodic solutions of the differential equation reach a limit when the second extremum is zero, since an arbitrary power of a negative functional value would necessitate a complex solution. Inspection of the series above indicates that the corresponding value for the first extremum occurs when each term is zero, or

$(c_0{)}_{\max }=(n+1{)}^{1/n}$

which is the limitation mentioned at the outset. This limit corresponds to an orbit with zero angular momentum.

The bracketing values of this limit can be found by taking a logarithm and applying l’Hôpital’s rule,

$\frac{\frac{d}{dn}}{\frac{d}{dn}}\frac{log(n+1)}{n}=\frac{1}{n+1}$

so that reversing the logarithm with an exponential gives

$\underset{n\to 0}{lim}(c_0{)}_{\max }=e\underset{n\to \infty }{lim}(c_0{)}_{\max }=1$

For values of the initial condition larger than this limit the constant l becomes negative, and being proportional to the square of angular momentum implies complex values for the latter. This potentially interesting feature of the equation will be left for a future presentation.

Formally integrating the differential equation provides an expression for the independent variable:

$t=\frac{1}{\sqrt{2}}\int \frac{dg}{\sqrt{g-\frac{g^{n+1}}{n+1}-l}}$

Integrating between the two extrema gives the time between, which as stated above is a quarter period of the Jacobi elliptic function. Replacing g with a dummy variable for clarity, the frequency of the Jacobi elliptic function in terms of the initial condition and the elliptic modulus is

$\omega =\frac{K\left(k\right)}{t_1}=\sqrt{2}K\left(k\right)[\underset{c_1}{\overset{c_0}{\int }}\frac{dc}{\sqrt{c-\frac{c^{n+1}}{n+1}-l}}{]}^{-1}$

The order of endpoints in the integral keeps the frequency positive. This is an exact expression for a given elliptic modulus k and ensures that the approximation always matches an exact solution of the differential equation at extrema. This can be seen by comparing the approximation so far with the numerical integration given at the outset:

Even though the denominator of the integral in the frequency is zero at both extrema, the integral can be numerically evaluated quite handily using tanh-sinh quadrature, apart from when c₀ = c₁ = 1 . In this limit integration is problematic and can be replaced with an exact result. If the extrema are written as c ≈ 1 ± δ , then applying the binomial series to the first two terms of the denominator gives

$\begin{array}{l}c-\frac{c^{n+1}}{n+1}\approx 1\pm \delta -\frac{1}{n+1}[1\pm (n+1)\delta +\frac{(n+1)n}{2}{\delta }^2]\\ \phantom{c-\frac{c^{n+1}}{n+1}}\approx \frac{n}{n+1}-\frac{n}{2}{\delta }^2=\frac{n}{n+1}-\frac{n}{2}(c-1{)}^2\end{array}$

Leaving the third term in the denominator initially untouched, the integral can be evaluated as an inverse sine. The constant l, being defined by the first two terms of the denominator at the extrema, is then replaced using the same expansion for the final result:

$\begin{array}{l}\underset{1-\delta }{\overset{1+\delta }{\int }}\frac{dc}{\sqrt{c-\frac{c^{n+1}}{n+1}-l}}\approx \underset{1-\delta }{\overset{1+\delta }{\int }}\frac{dc}{\sqrt{\frac{n}{n+1}-l-\frac{n}{2}(c-1{)}^2}}\\ =\underset{-\delta }{\overset{\delta }{\int }}\frac{du}{\sqrt{\frac{n}{n+1}-l-\frac{n}{2}{u}^2}}=\sqrt{\frac{2}{n}}{sin}^{-1}\left[u\sqrt{\frac{\frac{n}{2}}{\frac{n}{n+1}-l}}\right]{|}_{-\delta }^{\delta }\\ =\sqrt{\frac{2}{n}}{sin}^{-1}\left[\frac{u}{\delta }\right]{|}_{-\delta }^{\delta }=\sqrt{\frac{2}{n}}\pi \end{array}$

The approximation above in red is to the eye already a pretty good match for the numerical solution for arbitrary elliptic modulus. The discrepancy can be seen better by visualizing the difference between the numerical and approximate solutions:

Varying the elliptic modulus k shows that the approximation, while matching the numerical solution at extrema, is generally either slightly above or slightly below the latter between extrema. Careful variation with cursor keys will show that there is a small range of elliptic modulus for which the approximation appears to be equally above and below the numerical solution between extrema.

To make this more exact, define a simple deviation by

$\sigma =\sqrt{\sum _{\mathrm{evaluation\; points}}[g_{\mathrm{numerical}}\left(t\right)-g_{\mathrm{approximate}}\left(t\right){]}^2}$

where the summation is taken over some number of points equally spaced between extrema. The value of this deviation will of course depend the number of points used, but its value is less important than its general behavior. Choosing the number of points equal to ten for calculational expediency, deviation as a function of elliptic modulus looks like this:

Clearly there is a value of elliptic modulus for which the deviation assumes a minimum, at least for n ≥ 1 . The location of the minimum is insensitive to the number of points chosen.

A natural best choice for the elliptic modulus would thus appear to be the value that minimizes the deviation for a given exponent and initial condition. Estimating the behavior of these values as a function of the two other parameters will complete the approximate solution of the differential equation.

The behavior of the values is unclear for n < 1 because there is only an endpoint minimum at k = 0 . Since this single value is unlikely to give an accurate description of the solution over the whole range 0 < n < 1 , the remainder of the presentation will be restricted to n ≥ 1 .

Directly following this paragraph is a hidden interactive graphic ( show ) that evaluates locations of minima for given n as a function of c₀. Since the evaluations take some time, they have already been processed with this graphic and the resulting data embedded in the page. This allows subsequent visualizations of the data to be immediately interactive. The embedded data is explicitly constrained to be identically zero for both n = 1 and c₀ = 1 : the former condition matches the behavior of the exact solution for n = 1 , while the second is consistent with the behavior of the exact solution for n = 2 .

With the minima in deviation evaluated and embedded, visualizing their behavior is straightforward. In three dimensions this looks like

but is perhaps more usefully visualized in two dimensions as a function of exponent for a given initial condition,

or as a function of initial condition for a given exponent,

along with the behavior in exponent for the maximum allowed initial condition $c_0=(n+1{)}^{1/n}$ :

For most of their domains, the first two graphs are both reminiscent of an inverse secant, with the third mostly definitely so. All of these graphs go to zero for either n = 1 or c₀ = 1 , and since ${sec}^{-1}\left(1\right)=0$ the desired behavior is achieved by having an argument in (n - 1) and (c₀ - 1) offset by unity. The graphs also appear to have asymptotes to unity for large values, and this is achieved by dividing by the asymptotic value of the inverse secant. The behavior of the elliptic modulus that minimizes the deviation from the numerical solution is then of the form

$k(n,c_0)\approx \frac{2}{\pi }{sec}^{-1}[f_1(n-1)f_2(c_0-1)+1]$

where $f_1\left(0\right)=f_2\left(0\right)=0$ . Since the solution to this point is approximate and not exact, what is needed now are choices for these two functions that do a decent job of matching the available data. The functions should be as simple as possible, without worrying about finding a precise match to the data.

Begin with the third graph, the behavior in exponent for the maximum allowed initial condition, since it has fewer variables. One estimate arises from temporarily removing the dependence on c₀ by setting $f_2=1$ and letting the other function be a simple linear: $f_1=a(n-1)$ . Visualizing a variation of the single constant

indicates that there is a value between one and two where the two lines coincide quite nicely. A convenient way to set this value is from the exact solution for n = 2 , for which the limit value $c_0=\sqrt{3}$ and so

$\begin{array}{c}k=\frac{1}{\sqrt{2}}=\frac{2}{\pi }{sec}^{-1}(a+1)\\ a=sec\left(\frac{\pi }{2\sqrt{2}}\right)-1\approx 1.252\end{array}$

The remaining behavior for the argument of the inverse secant can be approximated by taking

$f_2=[\frac{c_0-1}{(n+1{)}^{1/n}-1}{]}^b$

The need for the exponent b can be seen by visualizing the dependence of elliptic modulus k as a function of exponent n along with the numerical data,

or as a function of initial condition for a given exponent along with the numerical data:

There again appears to be a value between one and two for which the estimate matches the numerical data fairly well, albeit not nearly as exactly as that for the behavior in exponent for the maximum allowed initial condition. Since there is no compelling way to select this value by comparison to known exact solutions, a value of b = 1.75 will be chosen. To the same level of approximation, let a = 1.25 .

With these approximate values for the parameters, the three-dimensional behavior in comparison with the numerical data is

which is a pretty good approximation to the surface, considering the complexity of the actual problem.

In summary, the differential equation

$\frac{d^{2}g}{d{t}^2}+g^n=1$

has a decent approximate solution of the form

$g\left(t\right)\approx c_1+(c_0-c_1){cn}^2(\omega t,k)$

with the three dependent parameters given by

$\begin{array}{l}{c}_1:\frac{c_1^{n+1}}{n+1}-c_1+l=0\\ k\approx \frac{2}{\pi }{sec}^{-1}[1.25(n-1)[\frac{c_0-1}{(n+1{)}^{1/n}-1}{]}^{1.75}+1]\\ \omega =\sqrt{2}K\left(k\right)[\underset{c_1}{\overset{c_0}{\int }}\frac{dc}{\sqrt{c-\frac{c^{n+1}}{n+1}-l}}{]}^{-1},l=c_0-\frac{c_0^{n+1}}{n+1}\end{array}$

In order for the solution to remain real, the initial condition lies between unity and

$(c_0{)}_{\max }=(n+1{)}^{1/n}$

Keep in mind that the behavior of the elliptic modulus was constructed to minimize the deviation of the approximation from the numerical solution. When applying this approximate solution to the orbits determined by power potentials with arbitrary real exponents, it will turn out that another choice is more appropriate.

Uploaded 2016.12.27 — Updated 2017.09.30 analyticphysics.com